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We present an efficient general approach to first principles molecular dynamics simulations based 
on extended Lagrangian Born-Oppenheimer molecular dynamics [A.M.N. Niklasson, Phys. Rev. 
Lett. 100, 123004 (2008)] in the limit of vanishing self-consistent field optimization. The reduction 
of the optimization requirement reduces the computational cost to a minimum, but without caus- 
ing any significant loss of accuracy or long-term energy drift. The optimization-free first principles 
molecular dynamics requires only one single diagonalization per time step and yields trajectories 
at the same level of accuracy as "exact" , fully converged, Born-Oppenheimer molecular dynamics 
simulations. The optimization-free limit of extended Lagrangian Born-Oppenheimer molecular dy- 
namics therefore represents an ideal starting point for a robust and efficient formulation of a new 
generation first principles quantum mechanical molecular dynamics simulation schemes. 



I. INTRODUCTION 

With the rapid growth of available processing power, 
first principles molecular dynamics simulations, where 
the forces acting on the atoms arc calculated on the fly 
using a quantum mechanical description of the electronic 
structure, are becoming an increasingly powerful tool in 
materials science, chemistry and biology pQ. While some 
early applications where performed already four decades 
ago [21 G] , it was not until the development of efficient 
plane-wave pseudopotential methods based on den- 
sity functional theory [TOl [H] and the fast Fourier trans- 
form [12] . that first principles molecular dynamics simu- 
lations became broadly applicable. 

There are two major approaches to first principles 
molecular dynamics: a) Born-Oppenheimer molecular 
dynamics [TH31 EHPJ and b) extended Lagrangian Car- 
Parrinello molecular dynamics [D1113131I31I1Q1I1QI1 
I17H19I . In Born-Oppenheimer molecular dynamics, the 
forces acting on the atoms are calculated at the relaxed 
electronic ground state in each time step, which provides 
a well defined and often very accurate approximation. A 
key problem, however, is that a straightforward imple- 
mentation of Born-Oppenheimer molecular dynamics is 
unstable and does not conserve energy without a high 
degree of convergence in the electronic structure calcu- 
lations. If this is not achieved, the electronic system 
behaves like a heat sink or source, gradually draining 
or adding energy to the atomic system H0J . Several 
techniques have therefore been developed that attempts 
to improve the efficiency of Born-Oppenheimer molecu- 
lar dynamics and reduce the computational cost of the 
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electronic optimization procedure [2T)H55] . In extended 
Lagrangian Car-Parrinello molecular dynamics, on the 
other hand, the computationally expensive ground state 
optimization is avoided. As in Ehrenfest based molec- 
ular dynamics [3SHSH], the electrons are instead treated 
as separate dynamical variables oscillating around the 
ground state. This approach permits a stable dynamics 
with a low computational cost per time step. Unfor- 
tunately, Car-Parrinello molecular dynamics simulations 
typically require shorter integration time steps and a 
system-dependent choice of electron mass parameters to 
yield reliable results in comparison on an "exact" Born- 
Oppnheimer molecular dynamics; although statistical av- 
erages are often in good agreement [TJ [TO] . 

Recently, an extended Lagrangian formulation for 
a time-reversible Born-Oppenheimer molecular dynam- 
ics was proposed [Ml [25], which combines some of 
the best features of Car-Parrinello and regular Born- 
Oppenheimer molecular dynamics, while avoiding some 
of their most serious shortcomings. It has been argued 
that extended Lagrangian Born-Oppenheimer molecu- 
lar dynamics can be seen as a general framework both 
for Born-Oppenheimer and Car-Parrinello molecular dy- 
namics |19) . In this modern formalism of extended 
Lagrangian first principles molecular dynamics, Car- 
Parrinello molecular dynamics appears in the limit of 
vanishing self-consistent field optimization |19j . How- 
ever, the optimization-free limit can be approached in 
different ways providing a variety of solutions. In 
this paper we show how extended Lagrangian Born- 
Oppenheimer molecular dynamics, in the limit of vanish- 
ing self-consistent field optimization, gives a first prin- 
ciples molecular dynamics at the same level of accuracy 
as "exact" Born-Oppenheimer molecular dynamics, but 
without requiring short integration time steps or a mate- 
rial dependent tuning of electron mass parameters as in 
Car-Parrinello molecular dynamics. The instability from 
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the systematic energy drift associated with incomplete 
convergence of the electronic structure in regular Born- 
Oppenhcimcr molecular dynamics is also avoided. Our 
work here represents a generalization and first principles 
extension of recent work that was demonstrated for semi- 
empirical self-consistent-charge tight-binding simulations 

The ability to achieve a high degree of accuracy in 
the limit of vanishing self-consistent field optimization 
serves two main purposes: 1) it simplifies the calculations 
with a reduction of the optimization cost to a minimum, 
and 2) it provides the ideal starting point for fully con- 
verged, i.e. "exact" , Born-Oppenheimer molecular dy- 
namics simulations when the requirement of accuracy is 
very high. The optimization-free limit of extended La- 
grangian Born-Oppenheimer molecular dynamics there- 
fore represents an efficient and robust framework for a 
new generation of first principles molecular dynamics 
simulations. 



II. EXTENDED LAGRANGIAN 
BORN-OPPENHEIMER MOLECULAR 
DYNAMICS 

Extended Lagrangian Born-Oppenheimer molecular 
dynamics [24) can be formulated in terms of a La- 
grangian, 

£ XBO (R, R, P , P ) = \ V MxR] - f/(R; D) 

1 (1) 
+ \^Tr[P 2 ]-^u 2 Tr[(D-P Q ) 2 }, 

where the regular Born-Oppenheimer Lagrangian defined 
at the electronic ground state density matrix D for a 
given set of nuclear coordinated, {Ri} = R, has been 
extended with auxiliary dynamical variables for the elec- 
tronic degrees of freedom, Po and Po, that evolve in a 
harmonic well centered around D. The potential en- 
ergy U (R; D) is here the Hartree-Fock or Kohn-Sham 
energy functional including the ion-ion repulsion energy 
|27j . The parameter /i is a fictitious electron mass and 
u! is the frequency determining the curvature of the har- 
monic well. Euler-Lagrange equations, in the limit /i — > 
[23] , gives the decoupled equations of motion: 



MiR! 



dU(R;D) 



dRr 



(2) 



Po =lu 2 (D~P q ). 



The partial derivative of U for the nuclear coordinate Ri 
is taken with respect to a constant Po, since Po is an 
independent dynamical variable. The equations of mo- 
tion can be integrated using a time-reversible symplectic 
scheme, both for the nuclear and electronic degrees of 
freedom [3SJ HO]. By using a time-reversible Pq as the 



initial guess of the iterative self-consistent field (SCF) 
optimization procedure, 



P, 



Pi 



Poo = SCF(Po) 



where 



D = lim D(P n ) = lim Z6 Ui Q I - Z T H(P n )Z) Z T 
the total Born-Oppenheimer energy, 



(3) 

1 

(4) 
(5) 



is stable without any long-term energy drift, even in the 
case of approximate convergence of P„ C|!HM1 [28 ■ The 
ground state density matrix D in Eq. Q is given from 
the Heaviside step function, 9, of the converged Fock- 
ian or Kohn-Sham Hamiltonian, i. e. for limn^,^ iJ(P„), 
in an orthogonal representation, Z T HZ, with the step 
formed at the chemical potential, fiQ , separating the occu- 
pied from the unoccupied states. The congruence trans- 
formation matrix Z is given from the inverse Cholesky 
or Lowdin factorization of the overlap matrix, S, deter- 
mined by ZSZ T = I . 



A. Fast quantum mechanical molecular dynamics 

As n — » in Eq. Q, i.e. in the limit of vanishing 
self-consistent field optimization, the equations of mo- 
tion for the extended Lagrangian formulation of Born- 
Oppenheimer molecular dynamics, Eq. ([2]), are given by 



M/P 7 



dU(TL;D(P )) 



OR! 



Po 



(6) 



Po = lo 2 (d(p )-p ) 



By avoiding the self-consistent-field optimization of Po, 
these equations of motion require only one single diag- 
onalization per time step in the construction of D(Pq) 
and therefore provide a computationally fast method for 
first principles quantum mechanical molecular dynam- 
ics (fast-QMMD) [55]. An alternative derivation of the 
fast dynamics represented by Eq. ([6| that is motivated 
through a different set of arguments is given in Ref. . 

To guarantee stability in the integration of the elec- 
tronic degrees of freedom in Eq. ^ , using an integration 
time step of St, the dimcnsionlcss integration parameter 
St 2 u> 2 , typically needs to be rescaled by a factor c € [0, 1] 
[26j compared to the original integration of extended 
Lagrangian Born-Oppenheimer molecular dynamics [30] . 
This stability condition further assumes convexity of the 
total energy functional between Po and D(Pq) [26] . 

The definition of D = lim^oo D(P n ) in Eq. Q and 
our particular choice of sequence of limits both for /i — > 
and n — > are important. For example, if we instead use 
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D = P n and let n — > in the Lagrangian (before deriv- 
ing the Euler-Lagrange equations of motion), we end up 
with a /x-dependent set of unconstrained Car-Parrincllo- 
like equations |19j and if ji — > already in the initial La- 
grangian, but with full self-consistency convergence, we 
recover (trivially) regular Born-Oppenheimer molecular 
dynamics. For our particular sequence of limits of /x and 
n, the fast-QMMD defined by Eq. ([6| is formally neither 
an extended Lagrangian nor a Born-Oppenheimer molec- 
ular dynamics. However, as will be demonstrated in our 
examples, the first principles fast-QMMD in Eq. ^ is 
a very close approximation of "exact" , fully converged, 
extended Lagrangian Born-Oppcnhcimcr molecular dy- 
namics. 



III. EXAMPLES 

A. Implementation 

Our fast-QMMD, Eq. ([6]), has been implemented based 
on Hartree-Fock theory in the Uppsala Quantum Chem- 
istry (UQuantChem) simulations package [55]) which is 
a freely available suite of programs for parallel ab ini- 
tio electronic structure calculations using Gaussian ba- 
sis sets, including Hartree-Fock and M0ller-Plesset per- 
turbation theory, configuration interaction, variational 
and diffusion Monte-Carlo, structural optimization, and 
first principles molecular dynamics. The nuclear coor- 
dinates are integrated using the velocity Verlet scheme 
and the electronic degrees of freedom with a modified 
Verlet algorithm, including a weak dissipation term to 
remove the accumulation of numerical noise |25[ I30j . 
Since P appears as a dynamical variable in Eq. a 
Hellmann-Feynman-like expression for the nuclear forces, 
under the constraint of Pq being constant, can be ap- 
plied. Thus, even if the ground state condition neces- 
sary for Hellmann-Feynman forces are not fulfilled, we 
still have a force expression of similar simplicity. For the 
basis-set dependent contribution we use the original ex- 
pression of the Pulay force term |31j . which provides a 
sufficiently accurate approximation |26j . Our first prin- 
ciples dynamics is implemented based on Hartree-Fock 
theory (27J E2] • The Hartree-Fock method is the start- 
ing point for correlated wavefunction methods and can 
be used as the computational prototype for density func- 
tional theory jTUl HU E3 EH] and hybrid schemes [3"S] . 
Our optimization-free Hartree-Fock molecular dynamics 
therefore demonstrates applicability for a broad class 
of first principles methods. Extensions to plane wave 
schemes should also be straightforward [25] , 



B. Molecular dynamics simulations 

Figure [T] shows the behavior of the total energy, Eq. 
([5]), for the simulation of a single water molecule using 
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FIG. 1: Total energy fluctuations for water using "exact" 
(5 SCF/step) Born-Oppenheimer molecular dynamics (XL- 
BOMD), Eq. pi, and the first principles fast-QMMD, i.e. 
XL-BOMD in the limit n -> 0, Eq. Q, in comparison to reg- 
ular Born-Oppenheimer molecular dynamics (BOMD), where 
the density matrix form the previous time step is used as the 
initial guess to the SCF optimization with the energy con- 
verged to < 0.01 ^iHartree. In (a) a STO-3G basis set was 
used, in the inset Eo = -74.949 au. In (b) a 6-31G" basis set 
was used, in the inset Eq = -76.0058 au 
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FIG. 2: Total energy fluctuations for a 130 ps simulation, us- 
ing "exact" (5 SCF/step) Born-Oppenheimer molecular dy- 
namics (XL-BOMD), Eq. g, and the fast-QMMD, Eq. {6}. 
Here E = -74.953 a.u. 
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FIG. 3: Total energy fluctuations for a 5 ps simulation of 
a small cluster containing 10 water molecules, using "exact" 
(5 SCF/step) Born-Oppenheimer molecular dynamics (XL- 
BOMD), Eq. Q, and the fast-QMMD, Eq. Here E = 
-749.774 a.u. The inset shows the temperature fluctuations 
for the first 200 fs of simulation. 




FIG. 4: Total energy fluctuations for ethane using "exact" 
(5 SCF/step) Born-Oppenheimer molecular dynamics (XL- 
BOMD), Eq. (2|, and the first principles fast-QMMD, i.e. 
XL-BOMD in the limit n -> 0, Eq. Q. In the upper panel 
(a) a time-step of At = 5 au was used. In the lower panel (b) 
a time-step of At = 10 au was used. Here -Bo = -79.224 a.u. 



a STO-3G basis set in (a) and a 6-31G** in (b). Reg- 
ular Born-Oppenheimer molecular dynamics, where the 
density matrix from the previous time step was used as 
the initial guess to the iterative ground-state optimiza- 
tion, exhibits an unphysical systematic drift in the to- 
tal energy because of the broken time-reversal symmetry 
[201 El]- This drift is avoided in the "exact" fully opti- 
mized extended Lagrangian Born-Oppenheimer molecu- 
lar dynamics (XL-BOMD), Eq. which is very close 
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FIG. 5: Total energy fluctuations calculated at a level of 
URHF theory for a H2O molecule, (a), and a CF4 molecule, 
(b), using "exact" (5 SCF/step) Born-Oppenheimer molecu- 
lar dynamics (XL-BOMD), Eq. pj, and the first principles 
fast-QMMD, i.e. XL-BOMD in the limit n -> 0, Eq. (6}. In 
the upper panel (a) a time-step of At = 20 a.u. was used. 
In the lower panel (b) a time-step of At = 40 a.u. was used. 
Here E = -76.043 a.u., in (a), and E = -429.57 a.u., in (b). 
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FIG. 6: Interatomic distances calculated at a level of URHF 
theory, using "exact" (5 SCF/step) Born-Oppenheimer 
molecular dynamics (XL-BOMD), Eq. pj), and the first prin- 
ciples fast-QMMD, i.e. XL-BOMD in the limit n -> 0, Eq. 

In (a), the interatomic distance between the Oxygen 
atom and one of the the Hydrogen atoms, Ko-h, in a H2O 
molecule. In (b), the interatomic distance between the Car- 
bon atom and one of the Fluorine atoms, Rc-f, in a CF4 
molecule. 
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FIG. 7: The root mean square deviation (RMSD) between the 
fast quantum mechanical molecular dynamics and "exact" (5 
SCFs/step) Born-Oppenheimer molecular dynamics, Eq. (J2|, 
for the nuclear forces (red squares) and for the total energy 
(green circles) calculated for four different molecules at dif- 
ferent time steps. For comparison the local error of the total 
energy (black filled circles) has been calculated. Simulations 
were performed with a 6-31G** basis set, using URHF theory 
as implemented in the UQuantChem code |29| . 
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FIG. 8: The deviation as measured by the Frobenius norm be- 
tween the fast-QMMD density matrix, Po, and the "exact" (5 
SCF/step) Born-Oppenheimer molecular dynamics, Eq. pj, 
(XL-BOMD) density matrix, D(Ps), after perturbing the den- 
sity matrix, Po(t), at t = 500 fs, by resetting Po at t — 500 fs 
to the initial density matrix at t = 0.5 fs. Simulations were 
performed for a single water molecule at room temperature 
with a time step of 0.5 fs using URHF theory as implemented 
in the UQuantChem code [29] . 



to the results from the optimization-free fast first princi- 
ples QMMD (red circles), Eq. as seen in the insets. 
In particular, any deviations between the optimization- 
free and the fully optimized Born-Oppenheimer molecu- 
lar dynamics simulations are small compared to the local 
truncation error, i.e. the amplitude of the total energy 
fluctuations that are caused by the finite size of the inte- 



gration time step St = 10 a. u. As in classical molecular 
dynamics, the dominating integration error is thus de- 
termined by our choice of integration scheme and the 
size of the time step. Figure [2] demonstrates the long- 
term stability of our first principles fast-QMMD, which 
shows no systematic drift in the energy over 120 ps of 
simulation time. However, for longer integration time 
steps we have occasionally noticed a small drift that 
seems to be caused by the dissipation force of the mod- 
ified Verlet scheme. This sensitivity, which not yet is 
completely understood, is not found in partially or fully 
SCF optimized versions of extended Lagrangian Born- 
Oppenheimer molecular dynamics. The next figure, Fig. 
|3j shows the corresponding simulation for a water clus- 
ter containing 10 water molecules simulated for a shorter 
simulation time. Because of the chaotic movements of 
the larger system, a direct comparison with respect to 
the total energy is harder. The inset shows a compari- 
son of the kinetic energy fluctuations given by the tem- 
perature over the first 200 fs of simulation time, shortly 
before they eventually get out of phase. The total en- 
ergy fluctuations of the fast-QMMD simulation shows a 
noisy behavior similar to a random walk compared to 
the "exact" Born-Oppenheimer simulation (XL-BOMD). 
Similar random walk-like noise have been seen in linear 
scaling XL-BOMD simulations [25] . 

Figure [4] shows the behavior of the total energy, Eq. (|5j, 
in simulations of a C2HQ molecule using a 6-31G** basis 
set, which represents a slightly larger and more complex 
system compared to the water molecule. As a comparison 
two different time steps were used, in panel (a) At = 5 
au and in panel (b) At = 10 au. 

In Figure [5j simulations of a H2O and a CF4 molecule, 
showing the behavior of the total energy, Eq. ([5]), using 
two times respectively four times as long time steps as 
the maximum time step used in the previous examples. 

Figurc[6]illustrates the behavior of the interatomic dis- 
tances, in simulations of a H 2 and a CF4 molecule. 

Figure [7] shows the convergence toward "exact" Born- 
Oppenheimer molecular dynamics as the length of the 
integration time step St is reduced. This scaling demon- 
strates how the fast-QMMD scheme provides a well de- 
fined approximation to exact Born-Oppenheimer molec- 
ular dynamics whith an error of order St 2 , i.e. 



dU(R;D(P )) 



OR! 



0{St 2 



(7) 



P =a J 2 (D(P )-P ) + O(St 2 ). 



The corresponding behavior was recently found in our 
studies based on self-consistent-charge tight-binding sim- 
ulations [26) . 

To illustrate the stability of first principles fast-QMMD 
we perturb a simulation by resetting the auxiliary density 
matrix Po(t) to its tg initial value after 500 fs of simula- 
tion time. During the continued simulation, the pertur- 
bation slowly disappears as Po(t) converges toward the 
electronic ground state, as seen in Figure [8j where the 
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deviation of the fast-QMMD density matrix Pq relative 
to the "exact" density matrix P5 is plotted as a function 
of time. This behavior demonstrates a key mechanism 
of our method. Instead of optimizing to ground state in 
each iteration as in regular Born-Oppenheimer molecular 
dynamics, the time evolution of the electronic degrees of 
freedom makes Po(t) converge toward the ground state 
dynamically. At convergence, the auxiliary density ma- 
trix Po(t) oscillates around the exact ground state with 
an amplitude that is of the order 8t 2 . 

IV. CONCLUSIONS AND SUMMARY 

The extended Lagrangian approach to first princi- 
ples molecular dynamics, as pioneered by Roberto Car 
and Michele Parrinello in its modern formulation of 
extended Lagrangian Born-Oppenheimer molecular dy- 
namics |191 124j , provides an efficient and versatile frame- 
work for first principles molecular dynamics simulations. 
Here we have shown how the ground state optimization 
requirement can be simplified and reduced to a mini- 
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